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ABSTRACT 

We examine the effects of density stratification on magnetohydrodynamic turbulence driven by the 
magnetorotational instabihty in local simulations that adopt the shearing box approximation. Our 
primary result is that, even in the absence of explicit dissipation, the addition of vertical gravity leads 
to convergence in the turbulent energy densities and stresses as the resolution increases, contrary to 
results for zero net flux, unstratified boxes. The ratio of total stress to midplane pressure has a mean 
of ~ 0.01, although there can be significant fluctuations on long (> 50 orbit) timescales. We find 
that the time averaged stresses are largely insensitive to both the radial or vertical aspect ratio of our 
simulation domain. For simulations with explicit dissipation, we find that stratification extends the 
range of Reynolds and magnetic Prandtl numbers for which turbulence is sustained. Confirming the 
results of previous studies, we find oscillations in the large scale toroidal field with periods of ~ 10 
orbits and describe the dynamo process that underlies these cycles. 
Subject headings: 



1. INTRODUCTION 

The magnetorotational instability (MRI) plays an im- 
portant role in determining the angular momentum 
transport rate (and therefore accretion rate) in most as- 
trophysical disks (Balbus & Hawley 1998). Therefore it 
is of considerable interest to understand what determines 
the saturation amplitude of the MRI in the nonlinear 
regime. Investigations of this question generally utilize 
numerical methods to study the time-dependent MHD in 
the local, shearing box approximation. 

Fro m the first three-dimensional studies (Hawley et al. 
Il995l ). it has been known that for uniform density 
the saturation amplitude depends on parameters such 
as the net flux and g eometry of the magnetic field 
threading the domain (|Sano et all [2004). More re- 
cently, there has been considerable interest in the ef- 
fect of microscopic dissipation, such as Ohmic resistiv- 
ity and Navier-Stokes viscosity, on the sa turation ampli- 
tude with various initial field geometries fe omang et al.' 
120071 : iLesur fc Longarettil [20071 : ISimon fc H awlev 2009). 
as well as the effect of the radial ex tent of the domain 
(Bodo et al. 20Qj; iJohansen erall [2009: Guan et al. 
[200I . 

One particularly important and puzzling result is 
that in the special case of no net magnetic flux 
with no explicit dissipation, the saturation ampli- 
tude of the MRI decreases w ith increasing resolution 
(jFromang fc Papaloizoul l2007l : iPessah et aTTl2007l ). In 
this case, it appears the amplitude of the microscopic 
diffusivities determines the saturation amplitude of the 
MRI. Although this result is of considerable interest from 
a theoretical perspective in understanding the MRI and 
MHD turbulence, it is not yet obvious it has application 
to real disks, in which the magnetic flux is unlikely to 
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be zero in every local patch for all time (Sorathia et al. 
2009), and which are vertically stratified. 

Although the saturation of the MRI has been studied in 
the local sh earing box approximation in vertically strat i- 
fied disks (jBrandenburg et all [T995[ : IStone et al.lll996f ). 
these early studies lacked sufficient computational re- 
sources to perform a systematic convergence study, or 
evolve the disk for hundreds of orbital times in order 
to measure accurately the saturation amplitude. In this 
paper, we use modern numerical methods to revisit the 
saturation of the MRI in vertically stratified disks^ with 
no initial net magnetic flux. Interestingly, in this case we 
find quite different results compared to the unstratified 
boxes studied by Fromang & Papaloizou (2007). In the 
stratified boxes studied here, the stress converges with 
numerical resolution even with no explicit dissipation. In 
fact, with explicit dissipation, we find in stratified disks 
there can be significant and sustained turbulence at mag- 
netic Reynolds nu mbers that suppress t he turbulence in 
unstratifed disks ([Fromang et al.l [20071 ). These results 
seem to be a consequence of an MHD dynamo that op- 
erates in stratified disks, and we explore the properties 
of this dynamo in this paper. 

This paper is organized as follows. In ^we summarize 
the most relevant properties of our numerical methods 
and describe our Fourier analysis. In ^ we report our 
results: the outcome of our resolution study in §3.1( the 
dependence on the vertical and radial aspect ratios in 
§3.21 an d §3.3( and the effects of adding finite dissipation 
in §3.4[ In 0we discuss the nature of the dynamo driving 
the sustained turbulence, and in ^ we summarize our 
conclusions. 



2. METHOD 

^ Throughout the paper we describe our simulations as stratified, 
even though we assume an isothermal equation of state. In this text 
stratification simply refers to the density stratification which is the 
result of vertical gravity in our equations. It is not a reference to 
the entropy gradient. 
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We use Athe na (|Gardiner fc Stand l2QQ5l . l2QQ8l : 
IStone et"al] l2QQ8l ) for all calculations presented in this 
work. We perform 3d MHD simulations, adopting the 
local shearing box formalism and including vertical grav- 
ity. We refer the reader to Stone & Gardiner (2009) for 
a detailed discussion of the equations, algorithms, and 
boundary conditions specific to the shearing box, as well 
as a description of their implementation in Athena. Here 
we just summarize the basic equations and the most rel- 
evant features for our current work. 

The local shearing box approximation adopts a frame 
of reference located at a fiducial radius corotating with 
the disk at orbital frequency Q. In this frame, the equa- 
tions of ideal MHD are written in a Cartesian coordinate 
system (x^y^z) that has unit vectors (i^j^k) as 



dp 
dt 



v-[H = o, 



(1) 



dpv 

~W 
dB 

~dt ' 



V- [pvv + J]=p^^(2qxi - zk) - 2Q.k x pv,(2) 



V X {v X B) = 0, 



where T is the total stress tensor 
T=(p+BV2)I 



(3) 



(4) 



I is the identity matrix, p is the gas density, p is the gas 
pressure, B is the magnetic field, v is the velocity and 
= B ' B. The shear parameter q is defined as 



<7 ^ 



dhir 



(5) 



SO that for Keplerian fiow q = 3/2. We adopt an isother- 
mal equation of state with p = c^p. 

In ^3. 41 we also present simulations which include terms 
for constant scalar viscosity and resistivity. The viscous 
term is V-M with 



pu 



dvi ^ dvj 2^..^ ^ 

dxj dxj 3 



(6) 



and the resistive term is — V x (77V x B) when added to 
the right-hand side of ([2]) and ([3|), respectively. 

These sets of equations admit the well know solution 
corresponding to (linearized) uniform orbital motion 



(7) 



which is used for the initial condition. The initial equi- 
librium density configuration is Gaussian with 



Po exp 



(8) 



where H = \/2cs/^ is the scale height of the disk. For 
consistency with previous work (Stone et al.l Il996l ). we 
take Vt = 10~^, = 5 x 10~^, and po = 1, yielding 
Po = 5 X 10~^ and H = 1. All simulations are initialized 
to have a weak magnetic field with a ratio of midplane gas 
pressure to magnetic pressure P = 2Pq/B'^ = 100. The 
configuration is a vertical field with zero net magnetic 
fiux that varies sinusoidally in the radial direction. 

We adopt boundary conditions which are shearing pe- 
riodic in X, and periodic in both y and z. Clearly, the 
periodic assumption in z is physically unrealistic in a 



stratified box. Of course, one advantage of this assump- 
tion is computational expediency, but the main motiva- 
tion is to give us some level of 'control' over the evo- 
lution of magnetic fiux in the simulation domain. Ver- 
tically periodic boundary conditions are useful in that 
the mean (volume averaged) toroidal field is conserved 
(i.e. remains zero). Outfiow boundary conditions will 
necessarily introduce electromotive forces (EMFs) at the 
vertical boundary which can drive growth of non-zero 
(By). Such mean field evolution is plausible in real disks, 
but we worry about spurious growth in {By) due to the 
manner in which outfiow boundary conditions are imple- 
mented. These considerations are important because the 
presence of mean azimuthal field may enhance or sustain 
turbulence (|Hawlev et al.l [19951 ). Since one of our pri- 
mary goals is to examine the robustness of MRI turbu- 
lence in stratified disks, this prescription, which prevents 
the grow of a (box integrated) mean field, represents a 
conservative approach. 

All simulations make use of Athena's orbital advec- 
tion scheme (Stone & Gardiner, 2009), allowing us to 
consider domain s with large radial extent. O rbital ad- 
vection schemes (|Masse tl l2000l : iJohnson et aII[2Q08 ) take 
advantage of the fact that Equations ([IE]) can be split 
into two systems, one of which corresponds to linear ad- 
vection operator with velocity Vk and another with only 
involves the fiuctuations dv = v — vk- The integration 
of linear advection operator is very simple and not sub- 
ject to a Courant-Friedrich-Lewy (CFL) condition. Since 
6v <C Vk near the boundaries, the CFL condition in the 
second system is much less restrictive than in standard 
algorithms, particularly for domains with large radial ex- 
tent. It also has the advantage of removing the system- 
atic variation of truncation error in troduced by the shear , 
which can lead to spurious effects (j Johnson et al.ll2008D . 

2.1. Fourier Analysis 

We utilize a number of diagnostic tools to analyze the 
simulation output, including Fourier analysis. This is 
straightforward in a periodic domain, but in a shearing 
periodic system, the basis functions are shearing waves 
with a time dependent wave vector. This complication 
can be handled with simple remap pings before and afte r 
Fourier transforming, as outlined in lHawlev et al.l (|l995l ). 

The quantities of principal interest will be the power 
density spectra (PSDs) of the magnetic and kinetic ener- 
gies. Although the PSD is highly anisotropic in /c-space, 
we still find it useful to plot shell averaged Fourier ampli- 
tudes. For example, we define the shell averaged power 
spectrum of the magnetic field as 

Bl^ink'imi^ (9) 
where \B{k)\'^ denotes the average of \B{k)\'^ over spher- 
ical shells, and B{k) = J B{x) exp {—ik • x)d^x is the 
Fourier transform oi B. ^ 

It is also instructive to look at the Fourier transform 
of the induction equation. Taking Fourier transforms of 
the X and y components of ([3]) we find 

+ (10) 



dt 



^ Of course, all Fourier analysis in this work refers to discrete 
Fourier transformations of discretized data. However, for the ease 
of readability, we will use continuous notation throughout the text. 
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TABLE 1 
Simulation Summary 



Simulation 


Domain^ 


Resolution 


Orbits 




{pVxSVy)/P0^ 


{By)/VP^- 


S32R1Z4 


H xAH xAH 


32/if 


300 


0.012 


0.0029 


0.066 


S64R1Z4 


H xAH xAH 


64/if 


300 


0.0075 


0.0018 


0.029 


S128R1Z4 


H x4H x4H 


128/if 


300 


0.0076 


0.0016 


0.034 


S32R1Z6 


H xAH x6H 


32/if 


165 


0.010 


0.0024 


0.074 


S32R4Z4 


AH xAH X AH 


32/if 


250 


0.0082 


0.0022 


0.035 


S64R4Z4 


AH xAH X AH 


64//f 


160 


0.0076 


0.0019 


0.040 



^ Lx X Ly X Lz 

^ Brackets denote temporal averages taken from 50 orbits onward and volume averages over whole 
domain. 

^ Brackets denote temporal averages taken from 50 orbits onward and volume averages over innermost 
two scale heights. 



and 



ld\By{k)l 



S -\- Ay -\- Ez^x + Ex^z- 



(11) 



2 dt 

We will focus on these two components as they appear 
to be the most important for understanding the disk dy- 
namo. 

The definitions of the terms on the right-hand sides of 
(fTO]) and ([II]) are 



S{k) = Re 
Mk) 
Ez^y{k)=Re 



Blik)- I Bx^^exp{-ik-x)d^x 



-Re 



dx 

Vsb^exp(-ik-x)d^x 
dy 



dy 



exp {—ik ■ x)d X 



(12) 
(13) 
,(14) 



and perform similar Fourier analysis on the three right- 
hand side terms individually, labeling them T^^, T55, and 
Tdiv, respectively. One drawback of this expansion is that 
terms such as B^Ovx/dx are present, even though they 
do not contribute to the evolution of because they 
appear with opposites signs in both Tdiv and Ti^y. Such 
terms can be quite large, complicating the interpretation 
of Tfo^, T55, and Tdiv We prefer to leave the right hand 
sides in terms of the EMFs. 

For plotting purposes, we find it useful to normalize 
the quantities on the right hand side of (p!Q|) and (pT|) 
with the power spectrum. To differentiate them from 
the unnormalized quantities, we will use lower case let- 
ters. For example, e^,^(A:) = 2Ey^z{k) / {\Bx{k)\^Vt). This 
then constitutes the Fourier amplitude of the normalized 
rate of field production of Bx due to the vertical varia- 
tion of £y The factor ^ has been introduced to make the 
quantities dimensionless rates. 



Ey^z{k) = -Re 

EzAk) = -Re 
and 

ExAk)=Re 



BE 

Bl{k) ■ I -^exp{-ik-x)d^x 



Kik) 



B;{k) 



dx 



exp {—ik • x)d^x 



exp (—ik • x)d^x 

dy 



(15) 



(16) 



,(17) 



where subscript i in Ai refers to either the x or y co- 
ordinate. The EMFs are defined as £^ = x B, with 
Vt = V — Vsh and 



LyL'^ 



II 



Vydydz, 



(18) 



where Ly and Lz the size of the computational domain 
in the y— and directions. The Ax and Ay terms are 
included for completeness, but they are generally much 
smaller than the other terms so we will not discuss them 
further. 

These relations are similar to the trans fer functions 
used i n iFroma ng ^ Papaloizou (2007) and ISimon et al.l 
([2009h . In fact,^ our definition of S is identical and if 
we sum Ai over all three spatial dimensions, it would 
equivalent to their definition of A. These authors expand 

V x{vtxB) = {B- \/)vt - {vt ■ \/)B - (V • Vt)B, 




Fig. 1. — Isosurface (at p = 0.75) and slices of the density at 
250 orbits in a domain of size AH x AH x AH (S64Z4R4). On the 
left face of the domain a slice of the magnitude of the magnetic 
field is shown. 



3. RESULTS 
3.1. Resolution Study 
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Our primary goal is to test the robustness of sustained 
turbulence and angular momentum transport in strati- 
fied shearing boxes with zero net flux, and in addition, to 
characterize the properties of the turbulence in this case. 
Figure 1 is an image showing the three-dimensional struc- 
ture of the density at late time (250 orbits) in a typical 
simulation, computed with a resolution of 64 grid zones 
per scale height in an AH x AH x AH domain. Spiral den- 
sity waves characteristic of all simulations of the MRI in 
shearing boxes are evident in the density isosur faces. 

0.05 on ' ' ' \ ' ' ' ' \ ' ' ' ' \ ' ' ' ' \ ' ' ' ' \ ' ' ' ^ 
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Fig. 2. — Sum of integrated Reynolds and Maxwell stresses as 
a function of time in if x 4if x AH stratified shearing boxes. The 
curves represent the 128/if (black, solid), 64/ H (red, dashed), and 
32/if (blue, dotted) resolutions. 

To investigate the convergence of the stress with nu- 
merical resolution, we have performed a resolution study 
at 32, 64 and 128 grid zones per scale height in an 
HxAHxAH domain (hereafter S32R1Z4, S64R1Z4, and 
S128R1Z4, respectively). Since interest in MRI turbu- 
lence is driven primarily by its role in angular momentum 
transport, we first focus on stress as a diagnostic. Table 
[1] and Figure [2] summarize the behavior of the stress as 
we vary the resolution. 

The MRI turbulence contributes to the stress through 
the Maxwell stress —BxBy and the Reynolds stress 
pVxSvy^ where Svy is the component of the velocity 
with background shear removed. The domain and time 
average of these quantities are listed in Table [H We 
have normalized them by the initial midplane pressure 
Pq. Since the initial condition is in hydrostatic equilib- 
rium and magnetic pressure is always significantly lower 
than gas pressure at the midplane (see e.g. Figure [6]), 
this value of the midplane pressure does not evolve signif- 
icantly. With this normaliz ation they are roughly eg uiy - 
alent to the a parameter of IShakura fc SunvaevI (|l973l ). 
The time average is carried out from 50-300 orbits to ex- 
clude the transient period of enhanced turbulence dur- 
ing and immediately after the linear growth phase of the 
MRI. 

Both contributions to the stress decrease as we increase 
resolution, but the changes are much smaller when going 
from 64/i^ to 128/i^ than from 32/i^ to 64/i^, indi- 
cating convergence. This should be compared wit h the 
beha v ior observed in unstrati f ied b o xes fe.g. Sano et al.l 
120041 : iFromang Papaloizoul l2007l : iPessah et all l2007f ) 



in which total stress decreases by factors of ~ 2 as 
resolution is increased by a factor of 2. The normal- 
ized total stress in the S128R1Z4 run is 0.0095, com- 
parable to prev ious results for stratified domains with 
zero net flux (Brandenburg et al."l995': 'Stone et al.l ll996l : 
..Johansen et a l. 2009; Suzuki & Inutsuka 2009). The 
Maxwell stress is 4-5 times greater than Reynolds stress, 
which is slightly higher than, but roughly consist ent with 
previ ous results for stratified domains (e.g. Sto ne et all 
Il996f ). Similar values are also observed in unstratified 
runs, where the result appears to be independent of field 
geometry or flux, an d depend mainly on the rate of shear 
(jPessah et alll2006f ). 

Table [T] also includes the rms toroidal field, vol- 
ume averaged over the central two scale heights and 
time averaged from 50 orbits onward. The rms field 
strengths are relatively weak, with (By)'^ ~ 0.01 (P^), 
but may still be dynamically important since the pres- 
ence mean toroidal field in unstratified simulations has 
been shown to enhance turb ulence stresses and energy 
densities (Ha wlev et al.l [19951 ). In fact, the rms toroidal 
field strength correlates well with the stress, although 
this may be the by-product of stronger turbulence rather 
than a cause. 
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Fig. 3. — Comparison of magnetic energy density power spectra 
for 32/H (solid), 64/if (dotted), and 128/H (dashed) resolutions 
in if X 4if X 4if stratified shearing boxes. Power spectra have been 
averaged over spherical shells of constant k = \k\ The normal- 
ization then makes the y ordinate proportional to the fractional 
contribution to the total power per logarithmic interval in k. The 
power spectra are time averaged from 50-300 orbits. 

Figure [2] shows the temporal variation of the total 
stress. There is considerable variability, with relatively 
long-lived (> 50 orbit) departures from the mean. In the 
S32R1Z4 run, the ~ 50 orbit periods of enhanced stress 
contribute almost as much to the average as the longer 
periods of 'quiescent' stress. There are similar long-term 
fluctuations in the higher resolution runs, but these are 
generally smaller in amplitude and less important for de- 
termining the overall mean. Nevertheless, it is clear that 
one must average over relatively long baselines to obtain 
a representative value, making higher resolutions pro- 
hibitively expensive. 

A power spectral analysis of the magnetic field also in- 
dicates convergence, at least for the large scales where 
most of the power resides. Figure [3] shows the averaged 
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PSD for the S32R1Z4, S64R1Z4, and S128R1Z4 simula- 
tions. We average over spherical shells in /c-space (see 
§2.ip and in time, excluding the first 50 orbits to avoid 
the initial transients. There is a significant drop in power 
in going from the to 64/ H resolution runs, com- 

bined with a shift in the peak of PSD to larger k. How- 
ever, when going from 64/ H to 128/ there is signif- 
icantly smaller drop in amplitude at most scales and a 
much smaller shift in the peak wave number, also indi- 
cating convergence. 
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Fig. 4. — Comparison of magnetic (solid) and kinetic (doted) 
energy density power spectra in the 128//f resolution, H x 4/f x 4if 
stratified shearing box. 



ated with the development, a nd subsequent buoy ant rise, 
of large scale magnetic fields (|Pessah et al.|[2QQ^ . 

In Figure [6] we show spacetime diagrams for several 
variables associated with the magnetic field. For the sake 
of brevity, we focus on magnetic quantities, which ap- 
pear to play the dominant role, as suggested by the PSD 
analysis above. The spacetime plots are generated by 
averaging over the x and y coordinates at each height in 
the domain every quarter of an orbit. The top panel of 
Figure [6] shows P = 2(P)/(5^), where the angle brackets 
denote horizontal averages. As noted previously, mag- 
netic pressure remains relatively weak near the midplane, 
but dominates in the surface regions {\z\ > 1.5H). 
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Although we focus on the magnetic energy density, a 
similar convergence is observed in the PSD of the kinetic 
energy density. We show a comparison of the kinetic 
and magnetic energy density PSDs for run S128R1Z4 
in Figure [H It is notable that power in the magnetic 
field fiuctuations exceeds that of the kinetic energy at 
all scales. This is in contrast to unstratified runs which 
typical show greater power in the kinetic energy at the 
lowest /c, with magnetic energy dominating the power 
at higher k. It also differs from simulations of helically 
driven turbulence where the kinetic and magnetic energy 
have comparable amplitude at all but the lowest k where 
magnetic energy dominates (Brandenburg "2001) . 

The convergent behavior of the stratified runs should 
be contrasted with that of the unstratified simulations 
shown in Figure [H This plot shows the PSD for 
four unstratified runs with resolutions of 32/i^, 64/ 
128/ H, and 256 /H in H x 4H x H shearing boxes. 
Each factor of two increase in resolution results in a 
decrease by nearly a factor of two in the integrated 
power. There is also a shift in the peak wavenumber 
to larger k as resolution increases. This lack of conver- 
gence is almost i dentical to the that found by previous 
authors (lF roman g"fc Papaloizoull2007l : iGuan et aLll2009l : 
ISimon etal.1 120091 ). It seems that both the power and 
characteristic scale of the turbulence are set by the do- 
main resolution. Adding stratification appears to fun- 
damentally change the dynamics and provides a charac- 
teristic scale and amplitude of the turbulence which is 
independent of the resolution. This could be related to 
the different mechanisms that lead to the saturation of 
the MRI in the presence of stratification, perhaps associ- 



FlG. 5. — Comparison of magnetic energy density power spec- 
tra for 32/if (solid), 64//f (dotted), 128//f (dashed), and 256/if 
(dot-dashed) resolutions in HxAH x H unstratified shearing boxes. 
Power spectra are time averaged from 60-100 orbits. The spectra 
peak at kH/{27T) = 4, 7, 14, and 17, for the 32/H, 64/if , 128/if, 
and 256 /H runs respectively. This roughly consistent with a scal- 
ing /cmax — 27r(if A)"-*^/^, where A is the spacing between grid 
zones. 

The middle panels show the normalized x and y compo- 
nents of the magnetic field. The By component is consid- 
erably larger than and they are negatively correlated. 
The symmetry of the boundary conditions and the initial 
conditions require the vertical average of these quantities 
to be zero. Nevertheless, there are localized regions of net 
field that, under the effects of buoyancy, trace out curved 
trajectories in the spacetime plot. This is similar to other 
'butterfiy diagrams' seen in pre vious shearing box calcu- 



1995; 



2009; 



lations of s t ratified domains [Brandenburg et alJ 

Stone et all Il996l: iTurnerl l2004l : iJohansen et aL. 

Suzuki fc Inutsukall2009l ). Consistent with previous sim- 
ulations, the polarity is usually even about the midplane 
and, at fixed height, oscillates quasiperiodically on time 
scales of < 10 orbits. 

Within the inner three scale heights, the horizontally 
averaged net field is a sum over a fiuctuating B field 
so that (By)'^ is much less than (By) and similarly for 
Bx. As these regions of net field buoyantly rise, the ratio 
of {By)'^/{By) increases. Near the vertical boundaries 

magnetic dominated regions of rather uniform B ~ Byj 
develop. Their presence is very likely related to our as- 
sumption of periodicity in the vertical boundary condi- 
tion, so they are likely not physically relevant to real 
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accretion flows. The degree to which these regions affect 
the dynamics is discussed further in §3.21 

The bottom panel of Figure [6] shows the Maxwell 
stress. In addition to the temporal fluctuations, there is 
also considerable variation with height. The middle three 
scale heights dominate and regions of greatest stress tend 
to be found off the midplane. The stress is weakest in the 
magnetically dominated regions very near the boundary, 
and is even slightly negative in some places (although 
the colorbar only goes to zero). A striking result is the 
correlation of regions of stronger than average {By) with 
regions of larger than average stress. This lends support 
to the idea that the presence of a mean toroidal field 
leads to enhanced turbulent stresses and energy densi- 
ties. Although not shown, we note that the spacetime 
plot of the Reynolds stress is very similar to the Maxwell 
stress and the two are well correlated in both z and time. 
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Fig. 6. — Spacetime plot for the 128/H resolution run in a 
H X AH X AH stratified shearing box. From top-to-bottom the 
panels show the horizontally averages of plasma /3, the normalized 
radial and toroidal components of the magnetic field (respectively), 
and the Maxwell stress as a function of height above the midplane. 



Figure [6] shows there is clearly a significant amount 
of large scale and long timescale structure in space and 
time coordinates, respectively. To better understand 
and quantify this structure, we perform a complimentary 
Fourier analysis in both space and time. Since we are 
interested in the structure of the horizontally averaged 



box, we focus on vertical /c- vectors with k = kzZ. Every 
one-quarter of an orbit, we compute the discrete Fourier 
transform B'^ikz) = A7:k1\B(kz)\^ ^ which is analogous to 
the B'^ defined in §2.H but with kz replacing k. Note 
that the two quantities can differ significantly since the 
Fourier amplitudes are far from isotropic in /c-space. For 
each kz^ we Fourier transform in time to obtain B'^(kz^ f) 
where / is the time domain frequency. We divide the 
data into five 50 orbit time series (between 50 and 300 
orbits), Fourier transform each separately, and then av- 
erage. Although we lose access to the longest timescales, 
the averaging reduces the variance in the resulting power 
spectra, which are plotted in Figure [71 




f (Orbit"^) 



Fig. 7. — Magnetic energy power spectra as function of both 
kz and /, where / is the frequency in the time domain. The top 
and bottom panels correspond to the amplitudes of the Bx and By 
contributions the magnetic energy density, respectively. A detailed 
description is provided in the text f N3.1|) . 



We plot both \og[fkzBl{kzJ)/{2Po)] (top panel) and 
\og[fkzB'^{kzJ)/{2Po)] (bottom panel). Note that the 
horizontal average (Bz) is conserved (at zero) to round- 
off error, so there is no physical information in the z 
component for vertical wave vectors. We have multiplied 
by both kz and / before taking the logarithm. Since we 
use a logarithmic scale for both the kz and / axes, this 
weights each pixel so that its color scales linearly with it 
contribution to the overall power (i.e. in the same sense 
that uFi, is commonly used in astrophysics). 

There are significant differences in the morphologies 
of the Bx and By power spectra. For large spatial scales 
(small kz)^ both Bx and By show a double peaked profile 
with significant power at large (~ 10 orbit) and small (< 
1 orbit) timescales, although the small scales are subject 
to aliasing. The dip at intermediate scales is somewhat 
more pronounced and persists to somewhat larger kz for 
By than for Bx. As we shift to larger /c^, the peak in 
Bx broadens significantly and the dip goes away entirely. 
For By there is a locus of maximum power which shifts 
to higher kz as / increases from about 0.1 to 1 cycles per 
orbit. 
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3.2. Vertically Extended Domains 

Due to the low densities and high magnetic field 
strengths near the vertical boundary that arise from 
stratification, increasing the vertical extent of the do- 
main is particularly computationally expensive. There- 
fore, we performed our resolutions study in boxes with 
four vertical scale heights, which seemed suitable to get a 
separation between the midplane dynamics and the ver- 
tical boundary. In order to confirm that our results are 
not strongly dependent on this choice, we have repeated 
our 32 /H resolution run with six vertical scale heights 
(hereafter S32R1Z6). As one can see in Table [TJ the 
resulting volume average of the stress in the two 
simulations is in reasonable agreement, although slightly 
smaller in the S32R1Z6 run. 
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Fig. 8. — Spacetime plot comparing the toroidal component of 
the magnetic field in stratified shearing boxes with 32/if resolution. 
The top and bottom panels show simulations with four and six scale 
height vertical extent (respectively). 

One downside to choosing vertically periodic bound- 
ary conditions, is the buildup of strongly magnetized re- 
gions with uniform By near the boundaries. Although 
these regions don't contribute significantly to the angu- 
lar momentum transport, they are likely unphysical so 
one might worry that they feedback on the dynamics 
closer to the midplane. In order to get a better sense of 
their effect on the fiow, we have plotted spacetime dia- 
grams comparing the S32R1Z4 and S32R1Z6 simulations 
in Figures ll] and [9l 

Figure [8] shows the y component of the magnetic field. 
In a larger domain, there are still regions of rather uni- 
form By near the vertical boundaries. They are some- 
what more extended in height but with a slightly weaker 
net field. The regions of net By generated near the mid- 
plane can buoyantly rise to larger heights before inter- 
acting with the boundary region. Since the horizontally 
averaged field tends to increase as the fluid rises, this lead 
to further enhancement of the field strength over those 
found in the smaller domain. 

Figure [9] shows the Maxwell stress in the two runs. 
In the central four scale heights the two plots look very 
similar, suggesting that the four scale height runs are 
yielding a fairly robust estimate for the angular momen- 
tum transport. Regions of enhanced Maxwell stress are 



again correlated with regions of strong net By. Similarly, 
the Maxwell stress is generally larger in the inner four 
scale heights than in the S32R1Z4 run. Note that the 
volume weighted average stresses in Table [T] are lower 
for S32R1Z6 than for S32R1Z4 because we are averag- 
ing over the whole box, including the regions of weak 
stress near the boundaries. If we restrict the averaging 
to the inner two scale heights for both the S32R1Z6 and 
S32R1Z4 runs, the time averaged Maxwell stress in the 
S32R1Z6 simulation is greater by about 10%. 
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Fig. 9. — Spacetime plot comparing the Maxwell stress in strat- 
ified shearing boxes with 32/if resolution. The top and bottom 
panels show simulations with four and six scale height vertical ex- 
tent (respectively). 



3.3. Radially Extended Domains 

For the sake of computational expediency, we have per- 
formed our resolution study on domains with only one 
scale height in the radial direction. This has tradition- 
ally been the box size employed in most shearing box 
computations, mostly due to CFL constraints on the 
timestep which are imposed by the background shear. 
Using Athena's orbital advection scheme (discussed in 
we can consider larger domains to examine the ef- 
fect of this choice on our results. We have computed 
shearing boxes with AH x AH x AH domains at 32 /H 
and 64/ H resolution (hereafter S32R4Z4 and S64R4Z4, 
respectively). 

We plot the evolution of the total stress in these two 
simulations in Figure [10] and the normalized, time and 
volume averaged stresses are listed in Table [H The mean 
values of the stress are very similar to each other and also 
to those found at higher resolutions runs in the smaller 
boxes (S64R1Z4 and S128R1Z4). This suggest that con- 
vergence is occurring at even lower resolution than in 
the smaller domain computations. The time evolution 
differs from that seen in Figure [2] in that the amplitude 
of fluctuations is much lower. 

In Figure [TT] we compare the PSDs from simulations 
with different radial extent but the same resolution 
(S64R1Z4 and S64R4Z4). The PSDs are very similar 
at all but the lowest k. At low k the differences are ac- 
counted for in part by our shell averaging scheme. Since 
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Fig. 10. — Sum of box integrated Reynolds and Maxwell stresses 
as a function of time in AH x 4H x AH stratified shearing boxes. The 
curves represent the 64/ H (black, solid) and 32/if (red, dashed) 
resolutions. 

the larger box is a AH cube, we can isotropically average 
shells all the way to k = 7t/{2H). Since we can not do 
this average isotropically with smaller box, some of this 
low k < 27r/H power is included in the k = 2t: / H bin. 
Overall, the PSDs seem to be rather independent of the 
aspect ratio, consistent with nearly identical values for 
the the box integrated stresses. 

Figure [12] shows the spacetime diagram for the 
S64R4Z4 run. Overah, it is rather similar to S128R1Z4 
spacetime diagram in Figure [6l There are long timescale 
(> 50 orbit) fluctuations in the Maxwell stress, as well as 
^10 orbit quasi-periodic variations in By and Maxwell 
stress which are qualitatively similar to those in Figure 
[6l However, the amplitude of fluctuations is generally 
smaller in the larger box, consistent with the lower am- 
plitude fluctuations in the total stress found in Figure 
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We flnd that the volume averaged stresses in each sub- 
domain are highly correlated with each other. Therefore, 
the lower amplitude is not simply the result of "averag- 
ing" over several independent boxes, but appears to be a 
global property of a correlated domain. This behavi or is 
somew hat surprising in light of the results of Guan e t al.l 
(|2QQ9l ). which show that the turbulence decorrelates on 
scales > in unstratifled boxes. 

To investigate the correlation in stratifled boxes, we 
have followed Gua n et al.l (|2009) and calculated the trace 
of the two-point correlation of the magnetic fleld 

iB = {SBi{x)SBi{x + Ax)) (19) 

where 5Bi = Bi — (Bi) and there is an implied sum- 
mation over i. We computed the average in two ways: 
using the full domain and using only the innermost two 
scale heights. The two different procedures give different 
results for the correlations in the x — y plane at large sep- 
arations, since full domain average is dominated by the 
magnetically dominated regions near the vertical bound- 
ary. Since these are likely unphysical, we report only the 
two scale height average. 
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Fig. 11. — Comparison of magnetic energy density power spectra 
in stratified shearing boxes with 64/ H resolution for domains with 
Lx = 4H (solid) and Lx = H (dotted). The larger radial extent in 
the former allows one to isotropically average shell to lower k. 

In order to better understand the lower fluctuation am- 
plitudes, we have split the larger domain into four sub- 
domains of one scale height each in the radial direction. 



Fig. 12. — Spacetime plot for the 64/ H resolution run in a 
4H X 4H X 4H stratified shearing box. From top-to-bottom the 
panels show the horizontally averages of plasma (3, the normalized 
radial and toroidal components of the magnetic field (respectively), 
and the Maxwell stress as a function of height above the midplane. 

We compute the correlation for 21 and 18 evenly 
spaced snapshots for for the S64R1Z4 (50-300 orbits) 
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and S64R4Z4 (50-150 orbits) simulations, respectively. 
Although correlation at small separations is relatively 
consistent from snapshot to snapshot, there is significant 
variation at larger scales. Therefore, we average over all 
snapshots to get a mean correlation for each run. We 
plot the results of the two scale height average for the 
X — y plane in Figure [131 For both simulations we nor- 
malize by their maximum values, ^o, which agree to 
within 10 %. Our results are qualitatively consistent with 
those of Guan et al] (|2009) in that we find comparable 
tilt angles 6t^ which is the angle between the major axis 
of the correlation and the azimuthal axis. The correla- 
tion is very similar in both simulations, although the tilt 
angle from the S64R4Z4 calculation is slightly smaller 
(l9t - 15° rather than 18°). 
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Fig. 13. — Two-point correlation functions for the magnetic field 
in the x — y plane. The axis labels x and y refer to the Ax and 
Ay implicit in (I19|) . We have normalized the by it's maximum 
value which occurs at Ax = Ay = 0. The left and right panels 
correspond to the S64R1Z4 and S64R4Z4 calculations, respectively. 

We also plot the correlation along the minor axis (de- 
fined by X cos Of -^ys'mOt) in Figure [Ml Again, the core 
of the correlation at small separations is nearly identi- 
cal in the two sim ulations and consistent with the re- 
sults of iGuan et al. (2009). Near Al - O.li^ the slope of 
the curve from the S64R1Z4 run fiattens and the large 
scale correlation plateaus at a value of — 0.04^o- 
The S64R4Z4 curve declines further, also fiattens with 
— 0.01^0 until Al ~ 1.3H where it drops nearly to 
zero. 

As equation (p!9|) requires, we have subtracted the 
mean field which is generally non-zero when only the 
inner two scale heights are considered (see Tabled]), but 
is identically zero when using the whole domain. These 
mean fields can be substantial and would dominate the 
large scale correlation if included. It seems that these 
fields are sufficient to enforce relative uniformity in the 
magnetic energy and Maxwell stress throughout the four 
scale height domain. This uniformity may disappear for 
sufficiently large domains, and we see some suggestion of 
this in a calculations at 32 /H resolution where we have 
compared SHxSHx 4H and 4Hx4Hx 4H domains. Al- 
though we find that in both cases the Maxwell stress and 



magnetic energy density in one scale height wide subdo- 
mains remain correlated, they show greater variance in 
the the eight subdomains of the large box than in the 
four subdomains of the smaller box. 
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Fig. 14. — Magnetic field correlation along the minor axis in the 
X - y plane, plotted for S64R4Z4 (solid) and S64R1Z4 (dotted). 
For the horizontal axis, Al is the displacement from Ax = Ay = 
along the minor axis of the correlation. 



3.4. Domains with Finite Dissipation 

It has been shown that the behavior of MRI turbu- 
lence in unstratified sh earing boxes depends on the na- 
ture o f the dissipation (ISano et al.l il998: Fleming e t al.l 
20001: iFromang et all [20071: iLesur Longaretti" 200^ 
bimon &: Hawlevll2009h . Simulations with explicit dis- 
sipation yield different results than those with only nu- 
merical dissipation, and the strength and evolution of 
the turbulence depends on both the viscosity u and 
resistivity r^, or equivalent ly the the Reynolds number 
Re = CgH/v and the magnetic Reynodls number Rm = 
CgH/r]. Specifically, it has been shown that Re and Rm, 
or alternatively the magnetic Prandtl number Pm = 
Rm/Re, determine whether turbulence is sustained in 
zero net fiux, unstratified s imulations ([Fromang et al.l 
[2007[lSimon Hawlevl[2009h . 

To see if these results hold in stratified domains, 
we perform three simulations with differing values of 
viscosity and resistivity: Re=800 with Pm=4 (here- 
after Re800Pm4), Re=800 with Pm=2 (Re800Pm2), 
and Re=1600 wit h Pm=2 (Rel6 0Pm2 ). Examination 
of Figure 11 of iFromang et alj (|2007f ) or Table 1 of 
iSimon Ha^^ J^OOS) 

show that none of these would 
sustain turbulence in an unstratified domains with no 
net field, regardless of whether Zeus or Athena is used 
for the simulations. We have confirmed these results for 
unstratified domains with our own Athena calculations. 

Figure [15] shows the stress as a function of time for the 
three simulations with a logarithmic vertical scale. The 
behavior of Re800Pm4 and Rel600Pm2 is rather differ- 
ent from the unstratified domains where the simulations 
decay rapidly to zero on timescales of 10-20 orbits after 
the initial linear growth of the MRI. Even Re800Pm2, 
which drops rapidly to a dimensionless stress of ~ 10~^ 
and then ~ 10~^ in the stratified domain, decays much 
more rapidly and continues to even lower values in the 
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unstratified domain. The behavior of the stratified do- 
mains is also considerably more complicated. Turbulence 
never decays away completely, but vigorous turbulence is 
not sustained in any of the calculations for longer than 
100 orbits. The amplitude of variability is large, and 
turbulence decays slowly on timescales of hundreds of 
orbits. The Rel600Pm2 and Re800Pm2 runs both show 
a recovery nearly to peak values after spending over 100 
orbital periods in stagnation or slow decay! 
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to a low amplitude, we have elected to exclude it. Of 
course, the amplitude of the power spectrum depends on 
the interval used in the time average, which is 25-100 or- 
bits. The Re800Pm4 and Rel600Pm2 power spectra are 
similar in shape, falling off somewhat more rapidly than 
S64R1Z4 as k increases. 

The power in Re800Pm4 in exceeds that in Rel600Pm2 
at all k. Note that these two calculations have the same 
resistivity but Re800Pm4 has a higher viscosity, so it 
has higher amplitude despite having larger overall dissi- 
pation (although this conclusion depends on the inter- 
val used for the comparison). Similar behavior is also 
observed in two-dimensional simulations of MRI driven 
turbulence with a vert ical magnetic fi e ld an d non-zero 
viscosity described in Masada fc Sand (|2008l ). In some 
of their runs saturation is not achieved since the tur- 
bulent stresses increase with time until the end of the 
runs. The fact that the magnetic field generated by 
the MRI can reach high amplitudes could be due to the 
viscous quenching of Kelvin-Helmholtz parasitic modes 
(iGoodman fc Xul TT994: Pessah & Goodman 2009). Al- 
though plausible, it is less obvious that this process is 
responsible for the similar behavior observed in the sim- 
ulations with stratification and non-mean magnetic fiux 
presented here. 



Fig. 15. — Sum of box integrated Reynolds and Maxwell stresses 
as a function of time in stratified shearing boxes with explicit dis- 
sipation. The curves represent computations with Re=800, Pm=4 
(black, solid); Re=1600, Pm=2 (blue, dotted) Re=800, Pm=2 
(red, dashed). 

Using the criteria of iFromang et al.l (|2007l ). we would 
probably have labeled the Re800Pm4 run as having sus- 
tained turbulence (over the first 100 orbits, which is the 
baseline used there), the Rel600Pm2 run as marginal, 
and the Re800Pm4 as either marginal or not having 
sustained turbulence, although the complex variability 
of the stratified runs makes this somewhat subjective. 
Figure 11 of Fromang & Papaloizou (2007) maps out a 
locus of sustained turbulence in Re - Pm space, and 
it's notable that Re800Pm4 and Rel600Pm2 simulations 
are on the cusp of showing sustained turbulence while 
Re800Pm2 more firmly in the non-turbulent regime. 
Therefore, it would seem that stratification slightly in- 
creases the parameter space for which sustained turbu- 
lence is possible, but does not qualitatively alter the con- 
clusion that turbulence dies out for sufficiently low Pm 
or sufficiently high Re. 

It is suggestive that all three sets of dissipation terms 
show sustained turbulence i n unstratiffed boxes on ce a 
net toroidal ffeld is imposed (|Simon Hawlevll2009l ). As 
noted previously, it is conceivable that main impact of 
stratiffcation is the production of toroidal ffeld, which 
then lead to enhanced turbulence. The turbulence in the 
stratiffed runs is signiffcantly less vigorous, but this may 
be consistent with the rms ffeld strengths in the stratiffed 
simulations b eing; much weaker than the toroidal ffelds 
considered bv lSimon fc HawlevI (|2009l ). 

In Figure [16] we plot the time average power spec- 
tra of the magnetic energy density from Re800Pm4 and 
Rel600Pm2, including S64R1Z4 for comparison. Since 
the magnetic energy in the Rel600Pm2 drops rapidly 
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Fig. 16. — Comparison of magnetic energy density power spectra 
for computations with (solid, dotted) and without (dashed) explicit 
dissipation. The explicit dissipation curves correspond to compu- 
tations with Re=1600, Pm=2 (solid) and Re=800, Pm=4 (dotted). 
Power spectra are time averaged from 25-100 orbits (with dissipa- 
tion) or from 50-300 orbits (without dissipation). 

4. DISCUSSION 

The numerical experiments presented here motivate 
two related questions: Why do the stratified simulations 
converge when the unstratified calculations clearly do 
not? What is the source of the dynamo cycles observed 
in the large scale fields? Answering the former requires 
detailed comparison with unstratified simulations, and 
will be the focus of future work. For the moment, we 
focus primarily on describing the large scale dynamo. 

First, we compare with the discussion of 
iBrandenburg et all (|l995l ). who noted the similarities of 
their stratified results with din a — ft dynamo model. 
They showed that the azimuthal EMF Sy = {Sv x B)y 
was correlated with the azimuthal field By in a manner 
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that leads to growth in the radial B^. Coupled to the 
shear, which regenerates By from Bx^ this describes a 
simple dynamo. 

Our results are in good agreement with their Equation 
(21). For the S128R1Z4 simulation, we find 

{£y) ^ (-3,5) X 10-^ {By) {Svy/^, (20) 

where the values in parentheses correspond to volume 
averages one scale height above and below the midplane, 
respectively. They also report a correlation of {£x) with 
{By) in their Equation (22). Again, our results are qual- 
itatively similar to theirs with 

{£,)c,{-l,l)xlO-^{By){5v^y/^ (21) 

above and below the midplane. The first correlation 
confirms that there is a mechanism for regenerating the 
poloidal field from a toroidal field. The second relation 
suggests that {£x) generally acts to reduce the magnitude 
of {By) at the midpl ane, and is (at least parti ally) the 
result of buoyancy, as Brandenburg et al.l (|l995 h discuss. 
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Fig. 17. — Power spectra (top) and normalized EMFs (bot- 
tom) for a single oscillation period of S128R1Z4. The top panel 
shows PSDs k\Ba:{ko)\'^/{2Po) (solid) and k\By{ko)\'^ /{2Po) (dot- 
ted), where the former has been multiplied by a factor of 400 for 
plotting convenience. The bottom panel shows ey^z(ko) (solid), 
ex,z(ko) (dotted), and s(ko) (dashed), which are defined in ^2.11 

It's instructive to examine this further with the Fourier 
analysis methods described in §2.11 We focus on the 
large scale field and consider the smallest vertically ori- 
ented vector k = k^z = 27r /LzZ. Since, kx = ky = 0, 
this term represents the Fourier amplitude of "horizon- 
tally averaged" quantities on the largest vertical scale. In 
Figure [T71 we plot time variation of magnetic fields and 
EMFs over a single dynamo cycle for this choice of k. 
The top panel shows the Fourier amplitudes of magnetic 
energy densities |5aj(/co)P (solid) and \By{ko)\'^ (dotted) 
for this wave vector. We have multiplied \Bx{ko)\'^ by 
a factor of 400 to plot both on the same scale. There 
is considerable variation from cycle to cycle, and Fig- 
ure [3 shows that these oscillations are only quasiperiodic 
with a broad range of power for periods near 10 orbits, 
depending somewhat on the wavenumber used for the 
analysis. Nevertheless, this example is typical in that 
the curves are out of phase with a more uniform varia- 
tion in \By{ko)\'^ than in \Bx{ko)\'^. 



The bottom panel shows the corresponding right hand 
side quantities in (p!Q|) and ([TTj), which, along with nu- 
merical dissipation, drive the evolution of the Fourier 
amplitudes. We normalize these quantities with the 
power spectra as described in ^2.1\ e.g. ey^z{ko) = 
2Ey^z{ko)/{\Bx{ko)\^n). We plot e^,^ (solid), e^,^ (dot- 
ted), and s (dashed). Note that ez,y{ko) and ez,x{ko) 
are zero for ko because of the periodic boundaries. The 
shear term s primarily drives the variation of \By{ko)\'^, 

fiipping sign as \Bx{ko)\'^ goes to zero. In contrast, Cx^z 
is generally negative, acting as turbulent resistivity. The 
ey,z{ko) term is more erratic, frequently flipping sign over 
a single cycle, but the net effect is an overall oscillation 
of \Bx{ko)\'^ over ~ 7 orbital periods. 

In many respects, the behavior we see in the stratified 
simulations is similar to that observed in the unstratified, 
zero- net fiux calculations of Lesur & Ogilvie (2008). Us- 
ing an incompressible spectral code, they find dynamo 
cycles with a ~ 5 orbit periodicity. This is similar to 
the oscillations in our stratified runs where rms power on 
large scales is broadly distributed on times scales ~ 6 — 10 
orbits (see Figure [7j). The normalized quantities plotted 
in Figure [iTl are equivalent^ to their Equation (16). Com- 
parison of Figure [T71 with Figs. 4 and 5 in their paper, 
show that the behavior of the EMFs during oscillations 
are also quite similar, suggesting that a common (or, at 
least, related) mechanism may be responsible for these 
oscillations. This motivates a more detailed comparison 
between unstratified and stratified runs in future work. 
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Fig. 18. — Time averaged and normalized EMFs for x (top 
panel) and y (bottom panel) components of the induction equa- 
tion. The curves correspond to the S32R1Z4 (black), S64R1Z4 
(blue), and S128R1Z4 (red) calculations. In the top panel we plot 
ez,y(k) (solid) and ey^z(k) (dotted). The curves in the bottom 
panel are ez,x(k) (solid), ex,z(k) (dotted), and s (dashed). With 
this normalization, curves for different simulations lie nearly on top 
of each other at low k. All quantities have been averaged from 50 
to 300 orbits. 

^ Their notation reverses the definition of x and y from that used 
here: their Bx and By are the toroidal and radial field components, 
respectively. 



12 



Although it is useful to focus on the vertical wave vec- 
tors when trying to understand properties of large scale 
fields, an understanding of the overall power spectrum 
benefits from an analysis of the shell integrated quanti- 
ties. We plot the time and shell average EMFs terms 
in Figure [m for S32R1Z4 (black) S64R1Z4 (blue), and 
S128R1Z4 (red). As in Figure [TTl these quantities are 
normalized by the shell integrated power spectra. Each 
normalized term is then time averaged from 50-300 or- 
bits. Since the amplitudes of the magnetic energy densi- 
ties seem to be in statistical steady states over this pe- 
riod, we presume the left hand sides of (p!Q|) and (pT]) are 
nearly zero. Therefore the sum of the terms in each panel 
must be balanced by num erical dissipation terms, as dis- 
cussed in pr evious work (jFromang fc PapaloizQuI l2007l : 
iSimon et al l [2009 ). 

In the top panel we plot the terms Cz^yik) (solid) and 
ey,z{k) (dotted) which contribute to evolution of B^. At 
the large scales, we find that the Cz^y term is more im- 
portant for field generation and its normalized amplitude 
is nearly independent of resolution. The Cy^z term is 
smaller in amplitude and slightly negative as large scales. 
Even though Cy^z tends to oscillates about zero over an 
individual dynamo cycle while Cz^y usually remains pos- 
itive, the amplitude of e^,^ is generally larger, so the 
dominance of Cz^y at large scales is not simply the result 
of time averaging. As one moves to smaller scales, Cy^z 
rises and eventually dominates the generation of Bx . The 
characteristic wavenumber at which the crossing occurs 
shifts to higher values as the resolution increases. 

The bottom panel shows ez,x{k) (solid), ex,z{k) (dot- 
ted) and s (dashed), the terms which contribute growth 
in By. At large scales growth of By is dominated by 
the shear term, while both Cz^x and Cx^z are of compa- 
rable magnitude and negative. At small scales, Cz^x and 
ex,z both grow, becoming positive and dominating over 
the shear term. Again, the wavenumber of the crossover 
increases with resolution. 

Authors have often focused on horizontally aver- 
age properties of the fiow or (equivalently) the power 
spec tral variation only along vertical wave vectors 
(e.g Fromang & Papaloizou||20j07|; [Lesur & Ogilvie""200^, 
which were discussed above). We note that the behav- 
ior of Cy^z and Cz^y we have described differs significantly 
from what one would infer if only vertical wavevectors 
were considered. As previously mentioned, the symme- 
tries of the periodic box force ez^y{kz) to be zero, and 
only ey^zjkz) contributes. However, it is clear from Fig- 
ure [181 that the vertical EMF and its toroidal variation 
is also essential for understanding the mechanism which 
sustains turbulence in these simulations. 

The question remains as to why the addition of stratifi- 
cation leads to convergence in the turbulent stresses and 
energy densities. One possibility is that development of 
local toroidal field is key to sustaining turbulence in both 
stratified and unstratified domains. It is possible that the 
strength of toroidal field is entirely set by the resolution 
in unstratified domains, while stratified domains offer a 
characteristic scale which is independent of resolution, 
due to the action of the large scale dynamo. Indeed, it 
has a lready been demonstrated (e.g. iHawlev et al .11 19951 : 
[Simon & Hawlevi 2009) that a global net toroidal field 
leads to enhanced turbulent energy densities and stresses, 
and leads to convergence in the stress as resolution in- 



creses (|Guan et al.ll2009l ). Furthermore, our simulations 
show a correlation between the stress and the strength 
of the mean toroidal field, both globally in the two scale 
height averages (Table [T]) and locally in the spacetime 
plots (Figures [6l and [T2|) . This hypothesis will be ad- 
dressed further in future research, comparing in detail 
the results presented here with those from unstratified 
runs both with and without mean fields. 

Due to our choice of periodic vertical boundaries, and 
our use of simplified thermodynamics, we have largely 
avoided detailed discussion of observational implications. 
Such questions are generally better addressed by studies 
which include more physically realistic v ertical bound- 
ary conditions (e.g. iMiller fc Stond [20001 ) or more real- 
istic thermodynamics, including the treatment of radia- 
tion (e.g. Turner 2004; Hirose et al. 2006). However, it is 
worth briefiy noting that our work co nfirms some impor- 
tant re sults of ear lier studie s fsee e.g. [Brandenburg et aP 
1995; Stone et al.l [l996: Mil ler Stonell2000[ ). Figures [3l 
and [6] show that a significant fraction of the magnetic 
energy in these simulations resides in large scale mag- 
netic fields that rise buoy antly to the low dens i ty reg ions 
above the disk midplane. iBlackman fc PessaE (|2009l ) ar- 
gue that the magnetic field structures that power accre- 
tion disk coronae must be associated with characteristic 
lengths that are large compared to the typical turbulent 
eddies. If this were not the case, the timescales associ- 
ated with turbulent diffusion would be smaller than the 
corresponding buoyant rise time, making it difficult to 
transport significant magnetic energy to the coronae. In 
other words, if the corona is a consequence of magnetic 
field structures that originate within the turbulent disk 
via the MRI (or other magnetic instabilities), but that 
dissipate above the disk midplane, then these structures 
must be of large enough scale to survive the buoyant 
rise without being shredded by the turbulence within the 
disk. Thus, the results presented in this paper provide 
support to the prevailing paradigm for X-ray emission in 
accreting systems which involves an optically thin, hot 
corona pow ered by the dissipation of magnetic fields (e.g. 
iHaardt fc Marasch i 1993; Field fc Rogers. 1993.). 

5. CONCLUSIONS 

We have used Athena to examine the effects 
of stratification on magnetohydrodynamic turbulence 
driven by the magnetorotational instability. We have 
shown that stratified simulations converge as reso- 
lution increases, even in domains with zero-net-fiux 
and no explicit dissipation. This is contrary to 
our own calculations of zero-net-fiux unstratified do- 
mains, which do not converge, confi rming previous re- 
sults (Froman gfc Papaloizoul 120071 : iGuan et all l2009l : 
ISh^ et al. 2001" We have also considered calculations 
with explicit dissipation, and confirmed previous results 
that the maintenance of sustained turbulence is mag- 
netic Prandtl number dependent. Stratification appears 
to extend the range for which sustained turbulence devel- 
ops, and may allow sustained turbulence at slightly lower 
Prandtl number for a given Reynolds number. However, 
the behavior is rather complex with larger variations and 
evolution on long timescales (greater than 100 orbits). 

At the highest resolutions considered {64/11 and 
128/i^) the ratio of total stress to midplane pressure has 
a mean value of a ~ 0.01, but with considerable fiuctu- 
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ation about this mean on long (> 50 orbit) timescales. 
Since real astrophysical systems are stratified, this some- 
what alleviates concerns that magnetorotational turbu- 
lence might be unable to provide the required angular 
momentum transport in accretion disks, although values 
a factor of ten hig her have been in ferred in some as- 
trophysical sources (|King et al.ll2QQ7f ). Similarly, it par- 
tially alleviates concerns that explicit dissipation may 
be required in global disk simulations at high resolution, 
as stratification and net toroidal fields arise naturally in 
such calculations. 

We have shown that these conclusions do not depend 
sensitively on the vertical or radial dimensions of the 
box7 Domains with radial extents of one and four scale 
heights give the same time averaged values for a and have 
nearly identical power spectral densities for the magnetic 
energy. Stresses are somewhat more sensitive to varia- 
tions in the vertical height of the domain, although this 
may be related to our assumptions of vertical periodic- 
ity. Increasing the vertical extent from four to six scale 
heights results in only a slight increase in the time and 
spatially averaged stresses, as long as the spatial average 
is carried out over the same volume (about 10% when 
using the inner two scale heights). 

Our results generally reproduce the qualitative features 
found by previous authors for stratified systems (e.g. 

^ Variations in the azimuthal length were not considered here. 



iBrandenburg et all Il995l : IStone et al.l I1996D . This in- 
cludes oscillations with a periods of < 10 orbits in which 
the horizontally averaged radial and toroidal fields alter- 
nate sign. Coupled with buoyancy this leads to a charac- 
teristic butterfiy diagram in horizontally averaged space- 
time plots. A compar ison of our results with those of 
iLesur fc Ogilvid (|2008l ) suggest the mechanisms for gen- 
erating the large scale field oscillations in the stratified 
and unstratified domains may be related. 
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